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ABSTRACT 

The fraction of star formation that results in bound star clusters is influenced by the density 
spectrum in which stars are formed and by the response of the stellar structure to gas expul- 
sion. We analyse hydrodynamical simulations of turbulent fragmentation in star-forming re- 
gions to assess the dynamical properties of the resulting population of stars and (sub)clusters. 
Stellar subclusters are identified using a minimum spanning tree algorithm. When consider- 
ing only the gravitational potential of the stars and ignoring the gas, we find that the identified 
subclusters are close to virial equilibrium (the typical virial ratio Qvir ~ 0.59, where virial 
equilibrium would be Qvii ~ 0.5). This virial state is a consequence of the low gas fractions 
within the subclusters, caused by the accretion of gas onto the stars and the accretion-induced 
shrinkage of the subclusters. Because the subclusters are gas-poor, up to a length scale of 
0.1-0.2 pc at the end of the simulation, they are only weakly affected by gas expulsion. The 
fraction of subclusters that reaches the high density required to evolve to a gas-poor state in- 
creases with the density of the star-forming region. We extend this argument to star cluster 
scales, and suggest that the absence of gas indicates that the early disruption of star clusters 
due to gas expulsion (infant mortality) plays a smaller role than anticipated, and is poten- 
tially restricted to star-forming regions with low ambient gas densities. We propose that in 
dense star-forming regions, the tidal shocking of young star clusters by the surrounding gas 
clouds could be responsible for the early disruption. This 'cruel cradle effect' would work in 
addition to disruption by gas expulsion. We suggest possible methods to quantify the relative 
contributions of both mechanisms. 

Key words: galaxies: star clusters - galaxies: stellar content - stellar dynamics - stars: for- 
mation - stars: kinematics - Galaxy: open clusters and associations: general 



1 INTRODUCTION 

Over the past years, the implications of clustered star forma- 
tion have touched a range of astrophysical disciplines, from 
the scales of the star forma tion process itself (see the review 
by iMcKee & Ostrike^ l2007h to the fundamental properties of 
young star clust ers (e.g. ' McMillan et al .l2007l:lAllison et al.l2009al : 
iMoeckel & Bon nell 200^, or po ssibly even the global stellar 
mass assembly of g alaxies (see e.g. lPflamm-Altenburg et alll2007l : 
ISastian et alj|2010f) . While it seems evident that mo st stars form 
in a clustered setting (e.g. IParker & Goodwinll2007l) . estimations 
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of the exact fraction are hampered by the substantial dissociation 
of stellar structure that occurs during (but is not necessarily re- 
lated to) the transition from the gas-embedded phase to classical, 
gas-p oor star clusters jLada & Ladal l2003l : IPortegies Zwart et al.l 
I2OI0I) . The traditional interpretation that most, if not all stars form 
in clusters, with gas ex pulsion leading to t heir early disruption 
('infant mortality', see |Lada & Lad3 20031 : ISastian & Goodwinl 
l2006l : lGoodwin & Bastianll2006t) has recently been challenged by 
observational studies suggesting that stars form with a continu- 
ous distribution of densities, of which only the high-density tail 
eventually leads to bound ste llar clusters l iBressert et alj i2O10l : 
iGieles & Portegies Zwarll201lh . 

Current advancements in numerical calculations of turbu- 
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lent fragmentation in star-forming regions enable the study of 
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cluste r ed star formation i n increasing deta i l (e. 
19981 : iKlessen & Burkerl) l200d : iBate et alj [2OO: 
20081) . However, theoretical investigations of the response of 
stellar structure to gas expulsion are still largely based on 
the assumption of either a static gas potentialj or initial 
equilibr ium between the stars and ^as (e.g. [Tu tukov 1978; 
'Adams 2000; Geyer & Burkert 2001'; Boily & Kroupi< I2003ij|bl; 
Baumgardt & Krou pa 2007; Parmentier et al. 200 8i) , which need 



not apply to star-forming regions in nature. A more realistic setting 
was recently explored bvl Offner et al.l ( 1200 9^, who find that the ve- 
locity dispersions of the stars in hydrodynamic simulations of star 
formation are smaller than that of the gas by about a factor of 5, 
suggesting that the assumption of equilibrium between both com- 
ponents indeed does not hold. The respo nse of star clusters to ga s 
expulsion has also been investigated by iMoeckel & Bate! ( l2010h . 
who consider N-body simulations of star cluste rs using initial con- 
dition s from hydrodynamic simulations, and bv lMoeckel & Clarkd 
JioT l|), who address the dynamical evolution of star clusters un- 
der the condition of ongoing gas accretion. They propose that the 
disruptive effect of gas expulsion is limited by the way in which 
gas and stars are redistributed by the accretion-induced shrinkage 

of clusters. 

The hydrodynamical calculations of iBonnell et al.l i2008t) 
cover the hierarchical formation of several st ellar (sub)clusters, 
which have been identified and analysed by iMaschberger et all 
( |2010() . The simulation is very suitable for investigating the proper- 
ties of the (sub)cluster population due to the relatively large range 
of mass scales of the modelled structure ( see Sect. O. In this pa- 
per, we analyse the simulations reported in lBonnell et alj l l2008h to 
probe the response of gas-embedded stellar structure to gas expul- 
sion. We consider the dynamical state of the stars while ignoring 
the gas, which is equivalent to observing the stellar structure un- 
der the condition of instantaneous gas expulsion at any time in the 
simulation. 

This paper starts with a discussion of the setup of the simu- 
lations, the subcluster identification algorithm, and the characteris- 
tics of the stellar structure in Sect.|2] An analysis of the dynamical 
state of the subclusters follows in Sect. [3] which covers quantities 
such as the virial ratio, bound mass fraction and gas content. The 
response of the clusters to gas expulsion and an extension to the 
length scales of actual star clusters is considered in Sect.|4] The pa- 
per is concluded with a summary and an outlook in Sect.|5] where 
we discuss the possible dependence of the results on the initial con- 
ditions of the simulations and the input physics, and suggest ways 
in which our analysis could be improved and extended. 



2 SIMULATIONS AND CLUSTER SELECTION 

In this work, we analys e the hydrodynamical/N-b ody simulations 
of iBonnelletal ' ('20(H) and Bonnell et al. limi), extending the 
analys is of IMa schberger et al. (2010) and Maschberger & Clarke 
( l201lh . These smoothed particle hydrodynamics (SPH) simulations 
follow the evolution of a initially marginally unbo und, homoge- 
neous gas sphere of lO' Mq with a diameter of 1 pc ( IBormell et al.l 
|200A . and a cylinder of 3 x 10 pc that contai ns 10"^ Ma gas, boun d 
in the upper part and unbound in the lower dBonnell et al.ll20()8h . 



' Except for a normalisation of the gas potential that decreases with time 
when the gas is expelled. 



Initial turbulent motions are modelled with an initially divergence- 
free random Gaussian velocity field with a power spectrum P(k) oc 
k''^ . The gas is kept at a temperature of 10 K, staying isothermal 
throughout the lO' Mq simulation. In the 10"* Mq s imulation, the 
gas follows a modified Larson-type equation of state jLarsonl2005h 
comprised of three barotropic equations of state: 



P oc p^, 

with P the pressure, p the density, and where 



7 = 0.75; P ^ Pi 

7 = 1.0; pi < p ^ p2 

7 = 1.4; P2 < p 5? P3 

7 = 1-0; p > p3. 



(1) 



(2) 



g cm 



5.5 X 10"'^ gcm-\p3 



and pi = 5.5 x 10" 
2 X 10^'^ g cm"'. For reference, a density of 1 Mq (0.01 pc)~' 
equals 1.6 x 10'" g cm"^ 

Star formation is modelled via sink particles, which are 
formed if the densest gas particle and its ~ 50 neighbours are 
gravitationally bound (the critical density for sink particle for- 



mation is 1.5 
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cm 
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for the 10 Mq si mulation 
Mq simulation, see lBormell et 



and 



6.8 X lO"'" gem' _ 

I2OO8I for details). The mass resolutions of the sink particles are 
~ 0.1 Mq and 0.0167 Mq, respectively. Accretion onto sink par- 
ticles occurs if SPH particles move within the sink radius (200 AU 
for both simulations) and are gravitationally bound, or if SPH parti- 
cles move within the accretion radius (40 AU for both simulations). 
Gravitational forces between sink particles are softened at smooth- 
ing lengths of 160 AU (lO' Mq) and 40 AU (10* Mq). 

Under the influence of gravity, the initially smooth gas dis- 
tributions quickly form filaments in which the sink particles are 
formed. The sink particles themselves group together in subclus- 
ters that merge into larger structures, leading to the formation of 
one 'star cluster' in the lO' Mq simulation and about three 'star 
clusters' in the 10"* Mq simulation. Throughout this paper, we will 
focus on the 10"* Mq simulation, which contains about ten times 
more subclusters than the lO' Mq simulation, and therefore en- 
ables us to consider a population of subclusters rather than a select 
set of examples. We also ran our analyses for the lO' Mq simu- 
lation, which gave results that are consistent with those from the 
10"* Mq simulation. 

For the identification of the subclusters, we employ a min- 
imum spanning tree (MST) based clustering technique, which 
has been used in the context of young star fo rming regions (e.g. 
IMaschberger et al.ll201Cl: [Kirk & Mversll201ll) . The MST, which 
has the advantageous property of not imposing any geometrical 
symmetry on the dat a set, has also been used to quantify the amount 
of sub s tructure (e.g.rCartwright & Whitworthll2004 ISchmeia et al.l 



20081: Masc hberger et al.l i201 0|) and mass segregation 



Allison et al. 2009a,b; Maschberger et alJ|2oTol: IParker et al.l 
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for alternative methods) in star forming regions. The MST is a con- 
cept from graph theory, which represents the unique connection of 
all points of a data set, so that there are no closed loops (a 'tree'), 
and so that the total length of all edges between the points is min- 
imal. Typically, two separated groups of points are connected with 
one single, long edge, whereas the points within the groups have 
short edges. By simply removing edges that are longer than a cho- 
sen break distance the tree can be split in subtrees, which connect 
the points of the subclusters in the data set (furthe r informatio n on 
MST based clustering techniques can be found in IZahnlll97ll) . To 
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Figure 1. Spatial distribution of sink particles that are present at the end of 
the Iff^ Mq simulation (t = 0.641 Myr), projected on the x-y plane. Black 
particles constitute subclusters, and dark grey particles belong to the field 
population. Since the spatial extent of the simulation in the z-direction is 
larger than in the x-y plane, some of the apparent clustering is the result of 
the projection. 



avoid spurious detections, we require that a subcluster contains a 
minimum number of 12 stars. 

The MST technique utilises one free parameter, which is the 
break distance. Automated methods to determine the break dis- 
tance are ill-suited for the analysis of the simulations due to the 
highly varying properties of the stellar distribution. We therefore 
choose a break distance of dbreak = 0.035 pc, which gives sub- 
clusters that are comparable to those identified by the human eye, 
although they do not include the stellar haloes surrounding them. 
This b reak distance is larger than the choice of iMaschberger et aU 
feOlOf) . but leads to comparable clusters because we analyse the 
stellar structure in three spatial dimensions and not in projection. 
Because the choice of a single break distance could introd uce an ar- 
tificial length scale into the analysis jBastian et alj|2007h . we have 
also performed our calculations for a set of other break distances 
in the range dbreak = 0.020-0.100 pc, of which the results are used 
when discussing the implications and applicability of our findings 
in Sects.|4]and[5] 

An example of the results of the subcluster identification 
method can be seen in Fig.[T] which shows the spatial distribution 
of sink particles and subclusters at the end of the lO"* M0 simula- 
tion at t = 0.641 Myr, some 0.3 Myr after the onset of star forma- 
tion. At this point in the simulation, after one free-fall time, about 
60% of all stellar mass is constituted by subclusters. The spatial 
distribution shown in Fig. [T] is the result of a complex tree of hi- 
erarchical merging of small subclusters ( Maschberger et al. 2010). 
This process is still ongoing at the end of the simulation, which 
is illustrated by the close proximity of the subclusters towards the 
right in the pla ne of Fig. [T] 

Following ICasertano & Huj ( Il985l) . we use the projected dis- 
tance to the Nth nearest neighbour to determine the surface den- 
sity distribution of the sink particles. For a rank number N, the 
local surface density at the locations of each of the sink particles 
is Hsink = (N — l)/(7rD^), with Djv the projected distance to the 
Nth nearest neighbour in the x-y plane. The resulting distribution 
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Figure 2. Surface density distributions of sink particles. The solid his- 
togram includes all sink particles present at the end of the simulation 
(t = 0.641 Myr), with the shaded histogram marking the subset of parti- 
cles that belong to subclusters. The dashed curve represents a lognormal 
fit to the distribution with a peak at log(5Isink) = 3.75 and a dispersion of 
""log 5: = 1.13. For comparison, the dotted histogram denotes the surface 
density distribution of sink particles at t = 0.442 Myr. 



of surface densities for N = 7 is shown in Fig. |2l which includes 
all sink particles at the end of the simulation (t = 0.641 Myr), as 
well as those from a snapshot at t = 0.442 Myr, not too long after 
the onset of star formation (also see Fig. [3](. The difference be- 
tween the surface density distributions at both times shows that the 
stellar structure in the simulation typically evolves towards higher 
densities, even though the density range spanned by the sink par- 
ticles does not change much. As should be expected, the high end 
of the surface density distribution is occupied by the sink parti- 
cles that are residing in subclusters (shaded area), reaching den- 
sities of more than 10^ stars per pc". These surface densities are 
several orders of magnitude higher t han those obse r ved in nearby 
(< 500 pc) star-forming regions by iBressert et alj j2010l) , which 
is not surprising for two reasons. Firstly, crowding obstructs the 
observation of the densest parts of star-forming regions, which are 
therefore not included in their sample. Secondly, the high densities 
that are achieved in the simulation are likely the result of the ini- 
tial conditions, with a mean initial gas mass surface density of over 
10'' Mq pc^^ in the x-y plane. Nonetheless, the surface densities do 
compare well to the high-density region in the Orion Nebula clus- 
ter (ONC), of which the central s urface density coincides with th e 
peak of the distribution in Fig.l2l( lHillenbrand & Hartmannlll998h . 
although the system under consideration here is about one order of 
magnitude younger than the ONC. 

The subcluster assembly history is considered in Fig. [3] which 
shows how the number of subclusters and the mean subcluster mass 
evolve as a function of time and the star formation efficiency (SFE, 
the fraction of the gas that has been co nverted to stars) of the en- 
tire simulation (see lBonnell et alfcoill for a detailed spatial analy- 
sis of the SFE). The number of subclusters initially increases until 
it reaches a maximum, which occurs when the formation of new 
concentrations of sink particles is neutralised by the hierarchical 
merging of other subclusters. This is nicely illustrated by the mean 
stellar subcluster mass, which steeply increases around t = 0.55- 
0.60 Myr, when the merging of new-formed subclusters causes the 
number of clusters to decrease. The mean mass keeps increasing 
until the end of the simulation due to the ongoing accretion of gas 
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Figure 3. Top: time evolution of the number of subclusters. Bottom: time 
evolution of the mean subcluster mass. The star formation efficiency (SFE) 
of the entire simulation is indicated along the top axis. 



and small subclusters onto sink particles and the merging of small 
subclusters with a few massive ones. This mass increase would 
eventually be halted on time scales much longer than the duration 
of the simulation, when the available gas reservoir is depleted or, 
more likely, when the further inflow of gas onto the subclusters is 
obstructed by feedback from supemovae and stellar winds (neither 
of which are included in the simulation). 



3 DYNAMICAL STATE OF STELLAR SUBCLUSTERS 

Whether or not a gas-embedded stellar structure survives gas ex- 
pulsion depends on its dynamical state. Excluding the gas from the 
dynamical analysis is equivalent to observing the stellar structure at 
the moment of instantaneous gas removal. This represents the ex- 
treme case of gas expulsion, since a more gradual expulsion could 
allow a subcluster to (adiabatically) respond to the potential change 
and thereby retain a larger number of stars. As a result, the analysis 
of the dynamical state of solely the stellar component in simula- 
tions of star formation provides a lower limit to the retention of 
stellar structure upon gas removal. 



radiufl the fraction of the subcluster mass that is bound, and the 
virial ratio. The properties of the stellar component are supple- 
mented with information on the gas, including the gas mass frac- 
tions within the subclusters. 

The gravitational boundedness and virial ratio of subclusters 
are fundamental measures for their dynamical state. Both quantities 
are based on the potential energy and internal kinetic energy of a 
subcluster. For a sink particle i, the potential energy V, is defined 
as 



V, = 



(3) 



where m, and m, are the sink particle masses and r,j their mutual 
distance. The kinetic energy T, of a sink particle is 



T, 



-m,\v, - Vd| , 



(4) 



where vi and Vci are the respective velocity vectors of the sink par- 
ticle and the centre of mass of the subcluster. A particle is gravita- 
tionally bound if T, -I- Vi < 0. We define the virial ratio (Jvir as 



(5) 



The factor 2 reflects the correction for counting the potential en- 
ergy twice for each particle pair when combining Eqs.[3]and[5] A 
subcluster is in virial equilibrium if Qvi, ~ 0.5 and gravitationally 
bound if Qvii < 1. Supervirial subclusters have > 0.5, while 
subvirial subclusters have (Jvir < 0.5. 

It is possible that a single, dynamically hard binary, triple or 
multiple system dominates the energy of a subcluster. We correct 
for this by searching the sink particle list for binarie|3and replacing 
them with a single centre of mass particle. We repeat this step two 
more times, thereby correcting for triples and higher-order multiple 
systems. During the last iteration, the kinetic and potential energies 
of the subclusters generally remain unchanged, which indicates that 
a correction for higher order multiples is not required. We quantify 
the effect of binaries on the observables of interest below. 

In previous studies, the SFE has been identified as a key pa- 
rameter which determine s the survival chances of s tellar clusters 
upon gas expulsion (e.g. iGoodwin & Bastiaij|20o3) . However, a 
more fundamental critical factor is the dynamical state of the stars 
when the gas is removed. The effective SFE, eSFE = 1 /2Qvii , was 
therefore introduced as a measure for the survival probability of 
stellar structure at the moment of instantaneous gas expulsion (e.g. 
[Verschueren 1990; Goodwin 2009). If the gas and stars were in 
virial equilibrium before gas expulsion, the eSFE is equivalent to 
the actual SFE. If they were not in virial equilibrium, the eSFE no 
longer reflects the actual SFE. In that case, the survival chance of 
stellar clusters is not related to the actual SFE, but is solely deter- 
mined by their dynamical state. The eSFE is naturally higher in the 
identified subclusters than in the simulation as a whole. 



3.1 Dynamical quantities 

We have followed the evolution of several (dynamical) properties 
of the identified sub clusters over the course of the simulation of 
iBonnell et al] i2008l) , such as the stellar mass, the stellar half-mass 



^ The subcluste rs have predominantly small elongations 

)Maschberger et alj |201(1 Fig. 10), which enables the use of a half- 
mass radius. 

^ Binaries are selected by identifying a most bound partner for each sink 
particle. If it exists and the semimajor axis is smaller than 1000 AU, it is 
considered a binary. 
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Figure 4. Time evolution of tlie mean virial ratio of tlie stellar .subclusters 
(black solid line), of the mean virial ratio weighted by subcluster mass (grey 
dash-dotted line), and of the total virial ratio of all sink particles in the 
simulation (dashed line). The horizontal dotted lines indicate the marginally 
gravitationally bound case (Qvir < 1) and the virialised case (Qvi, = 0.5). 
After t = 0.5 Myr, about 40-60% of the stars resides in subclusters for the 
break distance adopted here (d|,feak = 0.035 pc). 



Figure 5. Dependence of the bound mass fraction (ratio of the total mass 
of the bound sink particles to the subcluster mass) on the virial ratio for the 
subclusters from all snapshots of the simulation. Each symbol represents a 
subcluster. The dashed line represents a linear fit to the data. Like in Fig.O 
the dotted lines indicate the marginally gravitationally bound case (Qyi, < 
I) and the virialised case (Qyij = 0.5). 



3.2 Virial ratio 

As star formation progresses, the population of subclusters grows 
in terms of its total mass. The hierarcliical merging of the subclus- 
ters inhibits the increase of their number and causes it to level off 
towards the end of the simulation, when the formation of new sub- 
clusters is balanced by their accretio n onto more massive ones (see 
Fig. [3l and iMaschberger et"ai]|20 lOh . Another consequence of this 
hierarchical buildup is that the properties of the subcluster popu- 
lation as a whole is not a direct representation of the evolution of 
individual subclusters, but also include 'emergent' properties of the 
population due to the interactions between the subclusters. This is 
also relevant when considering the mean virial ratio of the sub- 
clusters as a function of time. Individual subclusters can be formed 
either subvirially or supervirially with respect to the total poten- 
tial, and would eventually virialise with the total potential if kept in 
isolation. Deviations from this trend occur when subclusters merge 
or accrete smaller stellar aggregates, which temporarily increases 
the virial ratio of the merger product due to the relative velocity of 
the progenitors. Another thing to keep in mind is that we only in- 
clude the stars in our dynamical analysis, implying that the obtained 
virial ratio is always higher than its actual value by an amount that 
depends on the gas fraction. This affects the mean virial ratio of the 
population of subclusters, in which there is a continuous formation 
of new, gas-rich subclusters, which are typically still supervirial 
and only reach Qvir = 0.5 after some further evolution. 

Despite the complex setting of hierarchical star formation, the 
evolution of the mean virial ratio can be used as a first indication of 
how the dynamical state of the subcluster population evolves over 
time. This is shown in Fig.|4l which also includes the time evolution 
of the virial ratio of the entire simulation. As indicated earlier, we 
ignore the contribution of the gas to the gravitational potential, in 
order to assess the dynamical state of the stellar structure under the 
assumption of instantaneous gas removal. Even without account- 
ing for the gas potential, the population of subclusters evolves to a 
near-virialised state on a time scale of only a few tenths of a Myr 
This would suggest that the subclusters are typically gas-poor on 
length scales corresponding to their half-mass radii. The difference 



between the mean virial ratio and the mean virial ratio weighted by 
subcluster mass in Fig. |4] indicates that more massive subclusters 
are typically somewhat closer to virial equilibrium than low-mass 
subclusters. This is more of a trend than a relation: a simple linear 
regression of the virial ratio and subcluster mass M yields a best fit 
of Qvir = 0.86 — 0.161ogM, but with scatter larger than the fitted 
slope. Lastly, it is also shown by Fig.|4]that the entire stellar popu- 
lation in the simulation does not reach virial equilibrium, but does 
become marginally bound. This occurs because the simulation as a 
whole has a higher gas fraction than the subclusters (see Sect. |3.3l >. 
which illustrates that the SFE depends on the location and length 
scale on which it is computed. The dynamical state of the entire 
simulation also bears some traces of the initial conditions, covering 
a cylinder that contains a bound upper half and an unbound lower 
half (see Sect.[2l(. 

The virial ratios of individual subclusters do not show no- 
table correlations with subcluster mass or half-mass radiufl In- 
stead, they depend more strongly on the recent mass evolution of 
the subclusters. The virial ratio temporarily increases whenever the 
subcluster mass increases, be it due to the merging with other sub- 
clusters or by individual sink particles moving inside the MST 
break distance. When sink particles move more than a break dis- 
tance away and the subcluster mass decreases, the virial ratio de- 
creases as well. Both are natural consequences of the inclusion or 
exclusion of transient substructure in the identification of the sub- 
clusters. The same dependence is found when using different MST 
break distances to identify the subclusters: larger break distances 
yield more extended subclusters and consequently the mean virial 
ratio is higher When set to extreme values (dbreak > 0.050 pc), 
close passages of subclusters are incorrectly picked up as merger 
products, causing a spurious increase of the virial ratio. For the 
largest break distance used in our analysis (dbreak = 0.100 pc), these 
fluctuations can yield mean virial ratios briefly hitting (Qvir) = 1, 
in clear contrast with the result from Fig.|4] 

The quantity that most strongly correlates with the virial ratio 

* Except for the unbound subclusters (Qvi, > 1), which are generally small 
(rh < 0.04 pc) and low-mass (M < 30 Mq). 
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Figure 6. Histogram of the virial ratios of the subclusters from all snapshots 
of the simulation (solid Une). The shaded histogram represents the set of 
subclusters from the last snapshot at t = 0.641 Myr. The dashed line is a 
Gaussian fit to the data for all snapshots, with mean value Qvir = 0.59 and 
standard deviation ctq = 0.16. The vertical dotted lines again indicate the 
marginally gravitationally bound case (Qvir < 1) and the virialised case 

(Qvir = 0.5). 



is the bound mass fraction of the subclusters, i.e. the fraction of 
their mass that is bound even without accounting for the potential 
of the gas (see the definition in Sect. l3.1l a sink particle is bound if 
Ti + Vi < 0.). It is shown in Fig.[5]that subclusters with high virial 
ratios tend to have lower bound mass fractions, albeit with substan- 
tial scatter. This is not surprising, because the virial ratio is effi- 
ciently increased by fast, unbound sink particles that are included 
by the MST algorithm but would be left out with a physically mo- 
tivated identification. For larger break distances, the correlation be- 
tween bound mass fraction and virial ratio is stronger, due to the er- 
roneous identification of kinematically hot structure as subclusters. 
Nonetheless, most subclusters contain only very few unbound stars, 
with typical bound mass fractions of 0.95. Figure|5]also confirms 
that most subclusters are close to virial equilibrium, which was al- 
ready suggested by the evolution of the mean virial ratio in Fig.|4] 
For further illustration, Fig. [6] shows the distribution of the virial 
ratios of the subclusters from all snapshots, as well as those from 
the last snapshot, which are shown as the shaded region. Only 8% 
of the subclusters are unbound when excluding the gas, while 25% 
remains subvirial. When considering the subclusters from all snap- 
shots, a Gaussian fit to the distribution of virial ratios gives a mean 
of Qvir = 0.59 and a standard deviation of (Jq = 0.16. As in Fig.|4] 
the gradual decrease of the mean virial ratio towards Qvir = 0.5 is 
also visible in Fig. |6] A comparison of the two histograms shows 
that the subclusters in the last snapshot are closer to being virialised 
than the population of subclusters from all snapshots. These virial 
ratios imply that the eSFE is close to unity, i.e. the majority of sub- 
clusters will not be strongly affected by gas expulsion (see Sect.|4ll. 



Figure 7. Histogram of the gas fraction within the stellar half-mass radius 
of each subcluster from the two snapshots at fi = 0.442 Myr and t2 = 
0.641 Myr (solid Hne). The shaded histogram only shows the gas fractions 
for the last snapshot at t = 0.641 Myr 



The bound mass fraction of the subclusters is not strongly 
affected by the presence of binaries (unbound sink particles are 
generally single), but because binaries are in virial equilibriun[f| or 
slightly subvirial, the mean virial ratio of the subclusters from all 
snapshots is decreased by 0. 1-0.2 if it is not corrected for multiples. 
About two thirds of this difference is due to binaries, while the re- 
maining third is accounted for by triples and quadruples. This shift 
of the virial ratio means that without correcting for binaries, the 
subclusters could be incorrectly interpreted as being slightly sub- 
viriaj3(Qvir ~ 0.45-0.50), and the entire simulation would be close 
to virialised (Qvir ~ 0.60-0.65) instead of the marginally bound 
state that is shown in Fig. |4] With respect to the binary-corrected 
results from Fig.|4l this rather modest difference arises because the 
finite gravitational smoothing length used in the simulation inhibits 
the formation of very hard binaries. Nonetheless, the correction for 
binaries improves the accuracy of our analysis, and therefore all re- 
sults shown in this paper are corrected for binaries and higher order 
multiple systems. 



3.3 Gas content 

The key question is why the subclusters are so close to virial 
equilibrium when neglecting the gas potential. An obvious answer 
would be that the subclusters are generally gas-poor, which would 
imply that they are hardly affected by the gas potential in the first 
place. To assess the gas potential and its time evolution, we have 
analysed the distribution of the gas in two snapshots of the simula- 
tion, at times ti = 0.442 Myr (when star formation is ongoing) and 
t2 = 0.641 Myr (the last snapshot of the simulation, after one free- 
fall time; also see Fig. [5). For each of the identified subclusters in 
these snapshots, we calculate the fraction of the total mass within 



Replacing hard binaries and higher order multiples by their 
centre of mass particles is essential to obtain a reliable picture of the 
subcluster dynamics. The disruption of the subclusters during gas 
expulsion is controlled by the dynamical state of the binary centres 
of mass rather than the binaries themselves. Had we not corrected 
for binaries or higher order multiples, the measures for the dynam- 
ical state of the subclusters would fluctuate with the orbital phase 
of a few tightly bound and eccentric binaries. 



' Instantaneously, this only holds for binaries on circular orbits. Binaries 
on eccentric orbits exhibit a vaiiation of the virial ratio, with a subviiial state 
near apocentre and a supervirial state near pericentre. When considering 
the time-averaged virial ratio, eccentric binaries are in virial equilibrium. 
However, because the phase velocity is lowest near apocentre, a sampled 
population of randomly oiiented binaries will typically be slightly subvirial. 
* As in all of Sect. 13.21 this statement excludes the gravitational potential 
of the gas. 
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the stellar half-mass radius of the stellar distribution that is con- 
stituted by gas. The distribution of these gas fractions is shown in 
Fig. [T] which confirms that the subclusters are indeed gas-poor on 
their typical length scales, with gas fractions of (fgas) = 0-0.2. Be- 
cause the simulation does not include feedback, this means that the 
accretion of gas onto the sink particles can keep up with the overall 
gas inflow towards the subclusters. Another mechanism that natu- 
rally leadsJo_gas;£oor subclusters i s their accretion-driven shrink- 
ing (iBomiell et aljl9 98: Moeck el & Batel2010l : lMoeckel & ClarQ 
which increases the density contrast between the subclusters 
and the surrounding gas. 

Gas accretion and the time evolution of the structural proper- 
ties of the population of subclusters both further decrease the gas 
fraction as time progresses. This evolution is illustrated by com- 
paring the data of the two snapshots in Fig. [T] corresponding to 
times ti = 0.442 Myr and t2 = 0.641 Myr. During the enclosed 
time interval, the mean gas fraction of all detected subclusters de- 
creases by 0.63 dex, from (fgas(ti)) = 0.238 to (fgUh)) = 0.056. 
The mean half-mass radius of the subclusterfij decreases from 
(ri,(ti)) = 0.020 pc to (rh(t2)) = 0.013 pc, which is a decrease 
of 0.18 dex. Even though the shrinking of subclusters is a second 
order effect caused by gas accretion, it is interesting to ask which of 
both mechanisms contributes most to the decrease of the gas frac- 
tion. Is it mainly driven by the increasing mean stellar density of the 
subclusters or by the ongoing gas accretion onto the sink particles? 

To assess the relative contributions to gas depletion by accre- 
tion and subcluster shrinking, we consider the spatial distribution 
of the sink particles and the gas. Due to the relatively small num- 
bers of stars in individual subclusters, it is best to examine the mean 
density profiles of the populations of subclusters in the two snap- 
shots at ti and t?. Such a combination of the different density pro- 
files decreases the influence of low-number statistics on the result. 
In Fig. [8] we show the subcluster mass-weighted, mean cumulative 
mass distributions of gas, sink particles, and both combined. The 
distributions represent the enclosed mass fractions /i, normalised 
to the sum of the subcluster mass Md and the enclosed gas mass 
within three stellar half-mass radii M„„s: 



Md -I- Msas 



(6) 



with i = {stars, gas, all} and Mi(^) the enclosed mass at ^ = r/ru, 
which is the radial distance in units of the stellar half-mass radius. 
The mean distributions shown in Fig.[8]are weighted by subcluster 
mass to emphasise those subclusters with better statistics. A first 
comparison of both panels in Fig. [8] shows that the gas fraction in- 
deed decreases between ti and t2. The contribution to this decrease 
by subcluster shrinking can be estimated by a simple thought exper- 
iment, in which the gas distribution is kept fixed and the distribution 
of stellar mass is compressed by the appropriate amount. In the top 
panel of Fig. [S] the half-mass radii {rh(ti)) and (rh(t2)) correspond 
to ^ = 1 and ^ = 0.66, between which the enclosed gas mass differs 
by 0.26 dex. In other words, if the gas distribution was held fixed 
and the stellar distribution was shrunk appropriately, then the gas 
fraction within the new half-mass radius would have declined by 
0.26 dex. This is a probe for the decrease of the gas fraction that 
is solely caused by the shrinking of the subclusters. Comparing it 
with the actual decrease of the mean gas fraction between ti and t2 
of 0.63 dex, we see that it covers about half of the decrease, with 



The subcluster sizes are strongly correlated with the adopted MST break 




0.1 



1.0 



distance. See Sect. l4. ll for a discussion. 



Figure 8. Subcluster mass-weighted, mean cumulative mass fractions 
(see Eq.|6) of sink particles (dashed line), gas pailicles (dash-dotted 
line) and the sum of both (solid line), as a function of the radial distance in 
units of the half mass radius (^ = r/ri,). Top: mean cumulative distributions 
for the subclusters present at ti = 0.442 Myr Bottom: mean cumulative dis- 
tributions for the subclusters present at t2 = 0.641 Myr. The shaded areas 
enclosed by the dotted lines mark the 16th and 84th percentiles and illus- 
trate the typical spread of the enclosed mass fractions of sink particles (dark 
grey) and gas (light grey). 



the remaining 0.37 dex covered by gas accretion itself - not only by 
adding to the mass in stars, but also by decreasing the gas mass. We 
conclude that the evacuation of the gas due to ongoing gas accretion 
is about equally important for the gas depletion as the shrinking of 
the subclusters. 

Apart from enabling a quantitative comparison of the effect 
of gas accretion and cluster shrinking. Fig. [8] also demonstrates 
the spatial variation of the gas fraction in the subclusters. At early 
times, the gas is still prevalent in the outskirts of the subclusters, 
contributing 20-60% of the enclosed mass at ^ = 3. At the end 
of the simulation this gas has mostly vanished, leaving only a few 
percent of the mass within the stellar half-mass radius as gas, and 
typically 20% at ^ = 3. The radial dependence of t he enclosed ga s 
fraction is qualitatively reminiscent of the model of lAdamj ( I2OOOI) . 
It is interesting to note that the relative increase of the enclosed 
gas mass fraction with respect to the enclosed sink particle mass 
fraction only occurs at radii where the latter flattens, i.e. the sub- 
clusters only become gas-rich at radii where very little stellar mass 
is present. The influence of the gas on the subcluster dynamics is 
therefore best evaluated at radii smaller than where the flattening 
of /isiars occurs. At ti , the ratio between the enclosed stellar mass 
and gas mass just before the flattening is about 4:1, while at t2 it 
has increased to 11:1. This suggests that if feedback starts at a time 
t > t2, the resulting gas expulsion will not strongly affect the sub- 
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cluster dynamics, and that their virialised state (see Sect. 13. 2t will 
be largely retained. 



4 RESPONSE TO GAS EXPULSION 

Motivated by the low gas fractions found in Sect.[3l we now address 
the response of the subclusters to gas expulsion in more detail. 



4.1 The expansion of subclusters 

The long-term response of the subclusters to gas expulsion can be 
evaluated by once again omitting the gas from the simulations and 
considering only the identified stellar subclusters and their evo- 
lution towards virial equilibrium. Given a certain virial ratio and 
bound mass fraction, does a subcluster expand or shrink after gas 
removal? We combine the data fr om the simulations with a simple 
energy argument similar to iHillsl jT980i) to estimate how the sub- 
cluster masses and half-mass radii evolve after the expulsion of the 
gas. It is insightful to consider the system at two key moments. 

(1) The time of instantaneous gas removal, which is equivalent to 
the current system in the simulations while omitting the gas. This 
can be done for each snapshot, thereby providing a larger sample of 
subclusters than when only the last snapshot were to be considered. 
Of course, including subclusters from different snapshots implies a 
correspondingly extended range of moments of gas expulsion. 

(2) The time at which each subchister attains virial equilibrium 
after removing the gas. These times are different for each cluster 
per definition, but by considering the subclusters at their respective 
times of virialisation the long-term impact of gas removal is most 
clearly isolated and identified. 

The evolution of the subclusters between these two moments can 
be quantified by evaluating the conservation of energy. For any sub- 
cluster, we can express the kinetic energy as 



T = iMc,(v-> = - VQ™- ^ ^Qvir, 



(7) 



where is the virial radius, and {v~) denotes the mean square ve 
locity, which as a result can be written as 

GMd 



(v-> = Qvi 



(8) 



The total energy at the moment of instantaneous gas removal thus 
becomes 



El = (Qvi 



1) 



2rv.i 



(9) 



where the relevant quantities have been marked with subscript ' 1 ' 
to indicate the moment of gas expulsion. 

Given a deviation from virial equilibrium, a subcluster will re- 
spond by changing its radius and/or mass. As can be verified from 
Fig.|5] most subclusters contain a certain number of unbound sink 
particles, which were either previously retained by the gas poten- 
tial, or are randomly passing the subcluster close enough to be in- 
cluded by the cluster identification algorithm. These unbound sink 
particles will escape the subcluster upon gas expulsion and take 
away some of the kinetic energy. We now consider a second mo- 
ment in time, at which the gas-rid subcluster has reached virial 
equilibrium (Qvir,2 = 0.5) and the unbound sink particles have suc- 
cessfully escaped. At this time, energy conservation dictates 



-E] = El + Eesc — - 



; +(Md,i — Md,2j — - — , 

4rv.2 2 



(10) 



where E2 is the total energy of the (virialised) subcluster, Eesc is the 
total energy of the escaped stars, and the relevant quantities have 
been marked with subscript '2' to indicate the moment of virialisa- 
tion. The parameter /? denotes the surplus energy per unit mass of 
the escaped stars after they clear the potential of the subcluster, in 
units of its mean square velocity. The values of /3 can be estimated 
from the simulation by computing (T, -1- Vi)/(GMd,2Ji3i/2rv,i) for 
each of the unbound sink particle^. For this, we use the re lation be- 
tween the virial and half-mass radius corresponding t o alPlummei 

l') potential, which is given by Tv = 1.3rh (e.g. .Fleggie & HutI 
2003). The escaping sink particles are the tail of an approximately 
Maxwellian velocity distribution of the sink particles in the sub- 
cluster, and consequently the distribution of /3 declines exponen- 
tially as f(/3) oc exp(— /3//3o), with /3o around unity. The mean of 
such a distribution equals /?o per definition, which illustrates that 
unbound stars typically retain velocities similar to the mean square 
velocity in the subcluster after they escape. 

Combining Eqs. [8]and[T0l one obtains an expression for the 
evolution of the gas-rid subcluster as it approaches virial equilib- 
rium, which relates the half-mass radii, masses, initial virial ratio 
and p. It is given by 



fh.2 
I'h.l 



J"v.l 



1 



1 - Qvir.I 



1+/3 



Md,2 

Md,i 



P Md.2 

2 Md.i 



(11) 



where rh,2 is the only unknown and all other variables given by the 
simulation. For /? = 0, i.e. all escaping stars are only marginally 
unbound, the expression returns the basic result that when unbound 
stars escape from a virialised system (Qvir.i = 0.5), it contracts to 
reattain virial equilibriurrQ. Inserting typical values of /3 = 1 (see 
above), Qvir.i = 0.59 and Md.2/Md,i = 0.95 (see Sect.[3]2ll in Eq.fTTI 
yields rh,2/rii,i = 1.04, which justifies the earlier assumption that the 
radius does not change much between instantaneous gas removal 
and virialisation (see footnote [8}. This minor expansion is driven 
by the slightly supervirial state of the subclusters, and inhibited by 
the escape of unbound stars, which have velocities larger than the 
escape velocity. 

As discussed in Sect. 13.21 the virial ratios of the subclusters 
give an indication of their survival fraction after gas removal. Out of 
all 140 subclusters identified in the simulation, only 10 have virial 
ratios Qvi, > 1 and are therefore unbound. In the last snapshot 
of the simulation, only one of the 21 subclusters will disperse af- 
ter the removal of the ga|3 As a result, typically 90-95% of all 
the identified subclusters survive gas expulsion. The fate of these 



° The denominator holds a slightly modified form of the mean squai'e ve- 
locity (v|) (see Eq.[8) and assumes that the virial radius does not change 
much between gas expulsion and virialisation (i.e. ry i ~ rv_2)- This is re- 
quired since ry 9 is not available in the simulation. The assumption will be 
verified below. 

' This situation, in which the naturally unbound component of a system 
escapes, should not be confused with the response of a viriaUsed system to 
mass loss due to stellar evolution, when energy is injected into the system 
to unbind mass, hi such a case, the surplus energy of the escaping mass 
is supplied by the energy injection and not by the dynamical system itself, 
which therefore mainly loses potential energy. This does not apply to the 
case under consideration in Eq. 1111 where no energy is injected and the 
unbound stars take away more kinetic energy than potential energy. 

These unbound subclusters ai'e typically low-mass, compact systems, 
and are often newly formed or have just experienced a subcluster merger. 
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Figure 9. Histogram of the stellar half-mass radii of the subclusters from 
all snapshots (solid line). The shaded histogram represents the set of sub- 
clusters from the last snapshot at t = 0.641 Myr Top: for the simulated 
radii, representing the moment of instantaneous gas removal. Bottom: for 
analytically computed radii, reflecting the subclusters at a later moment, 
when they have reached virial equilibrium. The dashed hnes are lognor- 
mal fits to the data for the subclusters from all snapshots, with mean 
values log(ri,/pc) = {—1.98,-1.95} and standard deviations criagi- = 
{0.73, 0.61} for the top and bottom panel, respectively, with median half- 
mass radii log(rii/pc) = {—1.96, —1.91}. The radii depend on the adopted 
break distance (see text), but the lognormal shape of the distribution per- 
sists. 

survivors depends on whether they expand, and how their environ- 
ment affects them. Expanded subclusters with lower densities are 
more susceptible to disruption by tidal shocks. The evolution of 
the half-mass radii after gas expulsion can be considered in more 
detail by evaluating Eq. [TT]for each of the subclusters in the sim- 
ulation. This enables a comparison of the distribution of half-mass 
radii of the current subclusters (at the moment of instantaneous gas 
removal) with the distribution of their half-mass radii when they 
have reached virial equilibrium, which is shown in Fig.[9l The dis- 
tribution of half-mass radii changes remarkably little after gas re- 
moval, as the means of the lognormal functions that are fitted to 
both distributions differ by 0.035 dex. This implies rh,2/'"h,i = 1-08, 
very similar to the earlier, simple estimate of rh,2/rh,i = 1-04. The 
subclusters in the last snapshot experience roughly 1.5 times this 
expansion after gas removal, which is of the same order of mag- 
nitude as the expansion of the other subclusters. As gas expulsion 
does not increase the cluster radii, they become no more suscepti- 
ble to tidal perturbations once the gas is removed than they are at 
birth. However, this does not exclude any expansion of the subclus- 
ters, since mass loss due to stellar evolution and tidal perturbations 
would still cause them to expand (also see Sect. 14. 3t . 

The characteristic size of the distribution in Fig. [9] strongly 
correlates with the break distance of the MST. As indicated in 



Sect. (2] we varied the break distance over the range dbicak = 0.020- 
0. 100 pc. By fitting lognormal functions to each of the obtained size 
distributions, we find that the corresponding mean sizes fh corre- 
late with break distance as fh/pc = (0.88 ± 0.16)(dbreak/pc) ' 
for 0.020 ^ dbreak/pc ^ 0.100. This demonstrates that the char- 
acteristic size of the size distribution has no physical meaning. 
Nonetheless, it is consistent with a lognormal distribution for all 
break distances, i.e. independently of the definition for the subclus- 
ter radius, and the smallest subcluster size does remain constant 
at about 0.005 pc, corresp onding to the smallest cores in the sim- 
ulation dSmith et alJliopgh . Regardless of the adopted break dis- 
tance, the smallest subclusters are always those with the lowest 
numbers of stars. These length scales are not very surprising, be- 
cause we are not considering actual stellar clusters, but the compact 
(sub)systems that develop during the early formation of a cluster. 

Indeed a length scale of a few 0.01 pc i s much smaller than 
that of typical embedded clusters (see e.g. iLada & Lada 200^ 
and of the order of the maximum size of Class protostars 
jLaunhardt et ai] |2010l) . The formation of multiple stars in such a 
small volume could be related to the use of sink particles in the sim- 
ulation, although it should be noted that small stellar clusters and 



multiple systems within star-forming 
on scales as small as 0.005-0.015 pc i 


globules are actually found 


Kraus & Hillenbrandl2008l; 


iGutermuth et alj2009l;lLaunhardt et al. 


201(](). Such separations are 



consistent with simulated low-JV subclusters of which the half-mass 
radii are dominated by a single massive star or binary. In Sect. 14.31 
the units of the simulation are rescaled to cover the typical mass, 
length and time scales of low-mass, young star clusters, causing the 
half-mass radii of the identified clusters to reach up to 0.3 pc. 

4.2 The cluster formation efficiency 

The instantaneous gas removal discussed in this paper is an ex- 
treme form of the more gradual expulsion occurring in nature. As 
a result, the described weak effect of gas expulsion should be even 
weaker in real subclusters. It seems that gas expulsion plays a neg- 
ligible role on the length scales of the compact stellar aggregates 
in star-forming regions. However, the regions between subclusters 
may still be gas-dominated, implying that feedback could prevent 
the further merging of subclusters and thereby inhibit their hierar- 
chical growth. 

The length scale on which the subclusters will have merged 
and have become gas-poor depends on the moment at which feed- 
back starts. For the MST break distance and corresponding length 
scale that is used in most of this paper, the subclusters are gas-poor 
irrespective of time. However, there should be a break distance at 
which a notable time-evolution of the gas fraction appears. By com- 
paring the subcluster gas fractions for different break distances, we 
find that at the end of the simulation (after one free-fall time), the 
subclusters have become gas-poor ((fgas) < 0.1) on a length scale 
of about 0.1-0.2 pc. This length scale will increase further with the 
number of free-fall times that are completed before the onset of gas 
removal. In turn, this increases the spatial extent over which the 
subclusters are allowed to merge before gas expulsion, which im- 
plies that the most massive bound structure is inversely related to 
the free-fall time. 

The free-fall time is related to the density as tff « p~^^^, 
which implies that the time of the onset of gas removal by feed- 
back tfb is associated with a density pfb that has a free-fall time 
equal to ta. For a g iven density spectrum of subclusters (see e.g. 
iBressert et al.ll201(]h . only those subclusters with densities p 2> 
Pfb oc have the opportunity to undergo the collapse and shrink- 
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age that we find in the simulations. The CFE increases with the 
fraction of subclusters that forms in these density peaks. As sub- 
clusters merge, accrete gas and shrink, the density of the stellar 
structure further increases (see Sect. O. Each free-fall time, more 
subclusters evolve into the density regime where p ^ pn,, also on 
larger length scales. This means that the length scales on which 
star-forming regions produce virialised stellar systems that are in- 
sensitive to gas expulsion are larger in dense sites of star formation 
than in sparse ones. The resulting dense clusters are also less sus- 
ceptible to disruptive tidal effects from their environment, which 
potentially further increases their survival chances. As a result, the 
CPE should increas e with density. Th rough the Schmidt-Kennicutt 
law jSchmidt|[T959l: lKennicutllll998h . this suggests a relation be- 
tween the CFE and the star formation rate per unit volume psFR 
or per unit surface area Tsfr- Indeed, first ob servational indica- 
tions for such a relation have been found by Larsen & Ric htleii 
( l200(]h .lLarsen (2004), and recently also by Goddard et al. (201^ 
who obtain CFE oc T^yr- ^ relation between the CFE and the 
star formation rate density would also be consistent with the high 
clust er formation effic iencies that are found in starburst galaxies 
(e.g. IZepfetaT]|l999 '). However, dense star-forming regions are 
generally also more disruptive to clustered s tructure due to the 
higher frequency and amplitude of tidal shocks jLamers et alj200^ : 
iKruiissen et al.l2011[) . This implies that a relation between the CFE 
and the ambient density would be weakened or could even be can- 
celled (also see Sect l4.3t . 

4.3 Infant mortality and the 'cruel cradle effect' 

It is often said that most stars form in st ellar clusters. The concept 
of 'infant mortality' jLada & Ladal2003l) . i.e. the rapid dispersal of 
stellar structure following the change of the gravitational potential 
due to gas expulsion, has been put forward in the literature to ex- 
plain the observed rapid, mass-independent decline of the number 
of stellar clusters between ages of a few Myr and several tens of 
Myr (e.g. lSastian & Goodwinll2006l) . Because the majority of stars 
has been thought to form in clusters, infant mortality was also held 
responsible for the low number of clusters per unit field star mass. 
However, recent (observational) evidence is pointing towards a pic- 
ture in which star clusters are the dense end of a continuous den- 
sity spectrum of star formation (see Fig.l2land lBressett et alj2oTol : 
iGieles & Portegies Zwartl201ll) . This view challenges the need for 
infant mortality in the early disruption of stellar structure. 

The results presented in this paper show that stellar substruc- 
ture can evolve towards a virialised state before the gas is removed. 
This occurs beca use the dynamics o f the stars and the gas are de- 
coupled (also see lOffner et al.l l2009^. as the accretion of gas onto 
the stars together with the subcluster shrinkage can compensate the 
overall gas inflow onto the subclusters. In time, this causes the sub- 
clusters to become gas-poor, thereby diminishing the disruptive ef- 
fect of gas expulsion. It depends on the length and mass scales on 
which the gas is evacuated whether this result can be extended from 
subclusters to actual star clusters. Towards the end of the simula- 
tion, after about 0.3 Myr of star formation, the subclusters have a 
mean mass of 40 Mq and are gas-poor on length scales of 0.1- 
0.2 pc. As a very crude first-order estimate, one can re-scale the 
units of the simulation to have the subclusters match the typical 
properties of young star clusterf^ By multiplying the mass unit by 
a factor of 25 and the length unit by a factor of 5, we re-scale the 

" The simulations are scale-free except for the sink particle radius and 



mean cluster mass to 10^ Mq, and the length scale on which the 
stellar structure will be gas-poor to 0.5-1 pc. By scaling the time 
unit accordingly, we see that such gas depletion is reached on a time 
scale of 0.8 Myr, which is of the same order as the expected tn, due 
to stellar winds and ionisation feedback. This order-of-magnitude 
estimate is of course far from conclusive, but it does show the rele- 
vance of pursuing this problem further. 

If clusters reach a relatively gas-poor state before the on- 
set of gas removal by feedback, the influence of gas expul- 
sion on the dynamical state of the clusters will be smaller 
than previously expected. Rather than leading to the disruption 
of clusters ('infant mortality'), the different spatial distributions 
of gas and stars imply that gas expulsion could leave clusters 
marginally affected, unbinding o nly a certain fraction of their 
stars (e.g. lMoeckel&Ba3l201(]|) . This is in clear contrast with 
earlier (theoretical) approaches in literature (e.g. Boilv & Krouia 
2003allbl: Isa stian & Goodwin 2006; Goodwin & Bastian 20o|; 
Baumgardt & Kroupj |2007| : iParmentier et all I2008D , which as- 
sumed a model where the gas and stars are in equilibrium dur- 
in g gas expulsion. Cle arly, this is not the case in the simulation 
of lBonnell eTai] | |2008|) . 

Because the density of the stellar structure determines whether 
or not gas expulsion affects the survival chances of star clus- 
ter s, a continuou s dens ity spectrum of young stellar structure as 
in iBressert et alj ( 1201 Ol ) and Fig. [2] naturally leads to the situa- 
tion in which the dispersed part of the new-born stellar structure 
is affected by gas expulsion, while the other, dense and clustered 
part is not. However, this does not imply that the survival chances 
of these clusters are necessarily higher. Recent papers have ar- 
gued that the disruption of star clusters due to tidal shocks from 
the n a tal environment could be substantial (Elmegr een & Huntej 
l20ld : iKruii ssen et all 1201 k Kruijssen & Bastian, in prep.). Al- 
though the disruption rate due to tidal shocks decreases with cluster 
density, sufficien tly strong shock s would still be able to disrupt 
dense clustersEI dCieles et alj|2006h . Such shocks could be preva- 
lent in dense star-forming regions. As clusters move out of their 

i irimordial environment, the typical ambient gas d ensity decreases 
Elmegreen & Hunteill2010l ; IKruiissen et al.ll201 ih . which lessens 
the disruptive effect of tidal shocks. Observationally, this mecha- 
nism affects the star cluster population in a way that is very similar 
to infant mortality: the fact that young clusters are more efficiently 
disrupted than older clusters gives rise to a strong decline of the 
number of clusters with age at young ages. This decline acts on the 
age scale corresponding to the time it takes to migrate out of the 
star-forming region. Rather than being an internal effect, like in- 
fant mortality is, the primordial disruption by tidal shocks is driven 
by the state of the environment in which the clusters are born. We 
will therefore refer to this mechanism as the 'cruel cradle effect'. 



accretion radius mentioned in Sect. |2] For the scaling used in this example, 
these respective radii are lO' AU and 200 AU. 

The strength of a tidal shock corresponds to the amount of energy it 
injects into the cluster 

Note that this early tidal di sruption of stellar clusters is very different 
from the mechanism studied bv lParmentier & Kroupj fcOllh . who consider 
the enhanced loss of stars due to the expansion of a young cluster to its 
tidal boundary that is induced by gas removal. The difference again lies in 
their assumption of equilibrium between the stars and the gas. If a young 
cluster has evacuated the surrounding gas and gas expulsion plays a minor 
role, then any gas expulsion-induced tidal stripping will be minor as well. 
Nonetheless, young clusters would still expand due to stellar evolutionary 
mass loss and the tidal shocks that are causing them to be disrupted. 
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It will be interesting to quantify what the relative contribu- 
tions of gas expulsion and the cruel cradle effect are to the early 
disruption of young stellar clusters. It is possible that both effects 
coexist, and that the relative importance changes with the environ- 
ment. It was explained in Sect. l4.2l that the fraction of clusters that 
is affected by gas expulsion decreases with the ambient density of 
the star-forming region. The cruel cradle effect shows the opposite 
dependence, as the disruptive effect of tidal shocks increases with 
the ambient gas density. It would therefore not be unlikely that gas 
expulsion and the cruel cradle effect each dominate a different side 
of the gas density spectrum of star-forming regions. Their relative 
strength would then determine the relation between the CFE and 
the ambient gas density. 



5 SUMMARY AND OUTLOOK 

In this paper, we have assessed the dynamical state of stellar struc- 
ture in star-forming regions and its response to gas expulsion by 
analysing the prope rties of the stellar structure i n the SPH/sink par- 
ticle simulations of iBonnell et al. Subclusters have 
been iden tified using a minimum s panning tree algorithm (MST, 
following iMaschberger et all l2010h . and binaries have been re- 
placed by their centre-of-mass particles when computing the global 
dynamical properties of the subclusters. We have also discussed the 
long-term implications of gas expulsion for the properties of star 
cluster populations. The main results are as follows. 

(1) The surface density distribution of sink particles follows an 
approximately lognormal distribution similar to that observed by 
[Sressert et al. (2010). However, the surface density corresponding 
to the peak of the distribution is several orders of magnitudes higher 
than the observed one, because the subclusters considered in our 
study are part of a region that would represent only one or two 
clusters in the observations. The high-density end of the distribu- 
tion is occupied by sink particles belonging to the subclusters that 
are identified with the MST. 

(2) When excluding the potential of the gas from the dynamical 
analysis and only considering the sink particles, we find that the 
simulation as a whole becomes marginally bound after one free- 
fall time, and the population of individual subclusters is close to 
virial equilibrium. The mean value of a Gaussian fit to the distribu- 
tion of virial ratios from all snapshots is Qvir = 0.59, where virial 
equilibrium would imply Qvir = 0.5. The mean virial ratio of the 
population slowly decreases with time, from Qvi, = 0.70-0.80 early 
on to Qvir = 0.55-0.60 towards the end of the simulation. 

(3) The virialisation of the subclusters occurs due to their low 
gas fractions. We consider the spatial distributions of gas and 
sink particles at two characteristic moments during the simulation 
(ti = 0.442 Myr and to = 0.641 Myr, reflecting the system early 
on and after one free-fall time), and find that the mean gas fraction 
within the stellar half-mass radii (rh ~ O.OI pc) of the subclus- 
ters decreases by 0.63 dex during the enclosed time interval, from 
{fgUh}) = 0.238 to (fg-Uh)) = 0.056. By comparing the density 
profiles of gas and sink particles, we conclude that this decrease is 
caused by gas accretion and subcluster shrinkage to approximately 
the same degree. 

(4) Because the subclusters are relatively gas-poor, they are only 
weakly affected by gas expulsion and the subsequent evolution to- 
wards virial equilibrium. According to our analytical estimate, they 
expand by an average factor of 1.08 after gas removal. The length 
scale on which the subclusters are gas-poor ((fgas) < 0.1) is about 
O.I pc at the end of the simulation. By scaling up the units of the 



simulation from subcluster to star cluster scales, we find that after 
about 0.8 Myr of star formation, star clusters with a mean mass of 
lO' M0 would be gas-poor on a length scale of 0.5-1 pc. 

(5) Only those (sub)clusters with densities much larger than the 
density that is associated with a free-fall time equal to the gas ex- 
pulsion time can exhibit the shrinkage and accretion that causes 
them to become gas poor. The fraction of clusters that reaches the 
required density to become insensitive to gas expulsion before the 
onset of gas removal therefore increases with ambient gas density. 
This suggests a relation between the cluster formation efficiency 
(CFE) and the ambient gas or star formation rate density, with a 
larger fraction of star formation resulting in bound clusters in dense 
regions. 

(6) A possible relation between the CFE and the ambient gas 
density is affected by a second mechanism. In dense regions, the 
survival chances of stellar structure are not determined by gas 
expulsion or 'infant mortality', but by the disruptive effect of 
tidal shocks from the surrounding gas felmegreen & Huntedl2OI0l : 
iKruiissen et alj[20I ih . This disruption of young clusters by their 
environment is titled the 'cruel cradle effect' and is suggested to 
take over the disruptive role of gas expulsion in dense star-forming 
regions. The strength and relative contributions of infant mortality 
and the cruel cradle effect as a function of ambient gas density will 
be the decisive factor to assess the relation between the CFE and 
the ambient gas density. This needs to be quantified in future stud- 
ies, because the time scale covered by the simulation is too short to 
include the disruption of subclusters due to the cruel cradle effect. 

Throughout the paper, we have mentioned directions in which 
further research should be undertaken to verify and expand our con- 
clusions. It is essential to check to what extent these results depend 
on the properties of the simulations we analysed, such as their ini- 
tial conditions and input physics. The three chief concerns would 
be whether the results are affected by (I) the turbulence spectrum 
and initial setup of the simulation, (2) the inclusion or exclusion of 
feedback and magnetic fields, (3) the choice of sink particle radius. 

(1) The turbulence spectrum and overall boundedness of the sim- 
ulation primarily influence the time evolution of the overall star for- 
mation efficiency ( McKee & Ostriker 2007; Dale & Bonnell 200^), 
and will only impact the evolution of subclusters if the inflow of gas 
becomes too high to be compensated by accretion and subcluster 
shrinkage. Judging the relative ease at which the subclusters in the 
current simulation become gas-poor, it is unlikely that the situation 
of a saturating gas inflow would take place. However, a dynamical 
analysis of a larger set of simulations will be needed to obtain a 
conclusive picture, also to include the formation of stars and sub- 
clusters in environments with lower densities. 

(2) Feedback from accreting sink particles would inhibit the in- 
flow of the gas, which in turn would lead to a lower gas fraction 
within the subclusters. There have been several efforts in literature 
to quantify the effect of (positive or negative) feedback on the star 
formation process, which generally consider effects on the length 
scales of the giant molecular clouds in wh ich the star formation 
takes place (see e.g. lMcKee & Ostrikej2007il While global effects 
could influence the gas-poor state of the (sub)clusters, the nature of 
feedback is such that it evacuates the gas on the spatial scales of 
the subclusters, which therefore should not lead to a fundamentally 
different conclusion than made in this paper. Magnetic fields could 
slow down the accretion process. If this decrease of the accretion 
rate would be smaller in the outskirts of the subclusters than within 
the stellar half-mass radius, it would increase the gas fractions and 
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virial ratios of the stellar component. Therefore, the influence of 
magnetic fields on our results will need to be investigated. 
(3) If the accretion and/or sink radii of the sink particles were 
comparable to the typical interstellar separation, the gas-poor na- 
ture of the subclusters would be a trivial result of a high 'filling 
factor' of the subclusters by the sink particles, as the volume where 
the gas could reside without being accreted would be too small to 
achieve a stable configuration. We have addressed this to first or- 
der by computing the accretion and sink volumes taken up by sink 
particles and dividing it by the enclosed volume. This was done for 
each sink particle while taking the nearest neighbouf^ and also by 
calculating a mean radial 'filling factor' profile for each subcluster, 
analogous to Fig. [8] Both methods returned low filling factors, with 
typical values of 10^" for the sink radius and lO^"* for the accre- 
tion radius. In other words, less than 1% of the volume inside the 
subclusters is taken up by the sink particles. In order to assess to 
which extent this allows us to neglect spurious accretion, we ran 
a set of simple test simulations with different accretion and sink 
particle radii. These tests show that the gas accretion rate is hardly 
affected by the accretion and sink radii, which validates our results. 
The details of the tests are given in the Appendix. 

Ideally, the next step would be to perform the same type of cal- 
culation as in Bonnell et al. (2008) for different initial conditions, 
including descriptions for radiative and mechanical feedback, po- 
tentially accounting for magnetic fields, and varying the sink parti- 
cle radius. The dynamical analysis of such simulations would pro- 
vide a good verification of our conclusions, and would improve the 
current understanding of the dependence on initial conditions and 
input physics. 

The order-of-magnitude extension of our results from subclus- 
ter to actual star cluster scales should be investigated further. With 
the continuously improving computational facilities, it will be pos- 
sible to simulate systems on the scales needed to cover the forma- 
tion of star clusters. The key ingredients of such an effort will be 
larger particle numbers and smaller sink particle radii. Addition- 
ally, infrared or spectroscopic observations can be used to verify 
the length scales on which star-forming regions are gas-poor prior 
to the onset of gas removal. The current and upcoming generation 
of telescopes will provide excellent opportunities for this. 

If gas expulsion indeed only weakly affects the survival 
chances of stellar structure, it will need to be verified in which 
regimes infant mortality still plays a role. In order to understand 
the relation between the CFE and the local environment, the rel- 
ative contributions to early star cluster disruption of infant mor- 
tality and the cruel cradle effect will need to be quantified. The 
kinematic signatures of both effects should differ and would there- 
fore be measura ble in the velocity maps of young disrupted clusters 
( lKruiissenl20TTl) . Possible ways in which this could be done obser- 
vationally include searching for young clusters that are currently 
going through gas expulsion and mapping the radial velocities of 
the stars, or tracing the velocity dispersion profiles of young, gas- 
poor clusters in dense regions. To aid this effort, the differences 
between the kinematic signatures of energy injection into a star 
cluster by gas expulsion or tidal shocks have to be established the- 
oretically. The combination of these approaches should provide a 
conclusive picture of the mechanisms that determine which frac- 
tion of star formation results in bound star clusters. 



We used the particle list that was corrected for multiples, implying that 
bound neighbours were ignored. 
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APPENDIX A: INDEPENDENCE OF RESULTS ON SINK 
PARAMETERS 

In this appendix we verify that the resolution of the SPH simulation 
is not playing an important role in the evolution of the stellar-to- 
gas mass ratio of the subclusters. We accomplish this via a series 
of controlled, idealised tests in which a cluster of 10 sink parti- 
cles accretes from an envelope of gas. The total mass of the system 
is 1, divided equally between the sinks and the gas. The sinks are 
initially of equal mass, thus each has mass 0.05. They are placed 
randomly in a Plummer model of virial radius fsinks = 1 , and we use 
the same initial configuration of the stars in each test. The median 
nearest neighbour separation of the sinks is 0.43. The gas is like- 
wise in a Plummer sphere spatially, although with a larger radius 
than the sinks. The gas has zero initial kinetic energy and minimal 
thermal support, so that the gas falls onto the sink cluster and is ac- 
creted. We run two sets of tests, one in which the gas sphere's virial 
radius is ten times rsinks, i.e. rgas = 10, and one in which rgas = 3. 

The two numerical scales we are concerned with are the 
accretion radius of the sinks face, and the smoothing length of 
the gas particles. For the sink radii, we use the set race = 
{0.125,0.0625,0.03125}. The middle value yields approximately 
the ratio of the neighbour distance to the accretion radius seen in 
the clusters in the simulation. The smoothing length of the gas is 
determined by the number of gas particles. To roughly match the 
simulated value, suppose the sinks have masses 1 Mq. The total 
gas mass is then 10 Mq, and 5 x 10'* gas particles approximates the 
resolution of the large-scale simulation. We run the 0.0625 accre- 
tion radius cases with four times more and fewer gas particles, i.e. 
2 X 10^ and 1.25 x 10*. 



Figure Al. Remaining gas mass as a function of time for our tests of the 
SPH simulation sink parameters. Top panel: the runs with rgas = 3rsinks; bot- 
tom panel, the runs with igas = 10rsi„iis. Grey lines have race = 0.125, black 
hues have race = 0.0625, and light grey lines have race = 0.03125. Solid 
lines have 5 X lO'' gas particles. The dashed lines have 2 X lO' particles, 
and the dotted lines have 1.25 X lO'' particles. 

In Fig. lAll we show the gas mass as a function of time for the 
test runs. In the top panel we show the results for the case with 
Tgas = 3. Time is measured dimensionlessly where we have taken 
the gravitational constant G = 1 ; the crossing time of the sink clus- 
ter is ~ 2. In this setup, the gas free-fall time is ~ 6 and the gas 
accretes quickly compared to the time for the N-body dynamics to 
dissolve the small-N sink system. The agreement between all the 
test runs is excellent. The bottom panel shows the rgas = 10 cases, 
and the gas falls onto the sink system over a longer time scale, with 
a free-fall time ~ 35. At early times the agreement is quite good, 
with some disagreement between the runs appearing after t ~ 20. 
We attribute this to the fact that at this point the N-body dynam- 
ics of the different runs have set the clusters on clearly divergent 
paths; recall that the gravitational smoothing length of the sinks is 
proportional to their sink radius. By t = 80 the cluster of sinks has 
effectively dissolved. We conclude that gas accretion time scale is 
not greatly affected by the choices of the sink radius. 

This paper has been typeset from a T^U ffl^ file prepared by the 
author. 



